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Abstract. The evolution of an isolated over-density represents a useful toy model to test the accuracy of a 
cosmological N-body code in the non linear regime as it is approximately equivalent to that of a truly isolated 
cloud of particles, with same density profile and velocity distribution, in a non expanding background. This is the 
case as long as the system size is smaller than the simulation box side, so that its interaction with the infinite 
copies can be neglected. In such a situation, the over-density rapidly undergoes to a global collapse forming a 
quasi stationary state in virial equilibrium. However, by evolving the system with a cosmological code (GADGET) 
for a sufficiently long time, a clear deviation from such quasi-equilibrium configuration is observed. This occurs in 
a time tLi that depends on the values of the simulation numerical parameters such as the softening length and the 
time-stepping accuracy, i.e. it is a numerical artifact related to the limited spatial and temporal resolutions. The 
analysis of the Layzer-Irvine cosmic energy equation, confirms that this deviation corresponds to an unphysical 
dynamical regime. By varying the numerical parameters of the simulation and the physical parameters of the 
system we show that the unphysical behaviour originates from badly integrated close scatterings of high velocity 
particles. We find that, while the structure may remain virialized in the unphysical regime, its density and velocity 
profiles are modified with respect to the quasi-equilibrium configurations, converging however to well defined 
shapes, the former characterised by a Navarro Frenk White-type behaviour. 
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1. Introduction 



In order to control the numerical integration accuracy of 
a finite self-gravitating system evolution in a static back- 
ground one may u se total energy and angular momen- 
tum conservation ( AarsethL 2003 ). In the cosmological 
case, instead of being conserved, the total (peculiar) en- 
ergy decreases with time. Indeed, the cosmic energy equa- 
tion, known as the Layzer-Irvine (LI) equation, describes 
the change of the total peculiar energy of a sys tem due 
to th e adiabatic expansion of the background ([Peebles, 
1983). The quantification of the deviations from the LI 
equation is one of the very first tests made in all cosmo- 
logica l codes (see, e.g., iTevssieil (|200ll) : iMiniati fc Colellal 
(|2007l) ). However, when typical initial conditions of stan- 
dard models of structure formation are considered, the 
application of the LI test is limited by several reasons. 
First of all, high redshift initial conditions of typical cos- 
mological N-body simulations are very uniform and cold, 
so that both the initial peculiar potential and kinetic en- 
ergies are close to zero: this makes difficult the numer- 
ical estimation of the cosmic energy equation. In addi- 



tion, at low redshifts, the density field is typically char- 
acterised by a great variety of non linear structures of 
different sizes and by low density regions still in the lin- 
ear regime. Structures of different sizes and in different 
physical regimes give very different contributions to total 
the peculiar potential and kinetic energies so that it is not 
straightforward to estimate which error in the cosmic en- 
ergy equation maybe tolerated and/or which constraints 
on the accur acy of the simulation maybe placed (see also 
discussion in Joyce &: Svlos Labinil ( 20131) and references 
therein). 

For this reason several authors have developed con- 
vergence tests to explore the role of the various physi- 
cal and numerical parameters that characterise a given 
simulation with the aim of establi s hing the reliability of 



the results (see e.g., iPower et al.l (j2003r) '). Alternatively 



there arc st udies that compar e differ e nt algorithms (see 
e.g.. [Knebe. Green fc Bennev I (|2003 ): lannuzzi fc DolasJ 
(|2011l ) and references therein). However, it is not a sim- 
ple issue to control the numerical accuracy of the non- 
linear gravitational clustering in a cosmological simula- 
tion; for instance the question of determining the optimal 
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spatial resolution, i.e. the best value for the gravitational 
force softening lengt h, is still an open pro b lem (see e.g., 
Splinter et alTdigSSl): IPower et all (l2003l): iRo meo et al 



(l2008[ ): lJovce. 
In this 



Marcos fc Baertschiger 
pape r, fol lowing 



2OO9I)). 

the work of 
consider the evo- 
namely an isolated 



Joyce fc Svlos Labinil (|2013l ). we 
lution of a very simple system, 
over-density in an expanding background, as a toy model 
to control the accuracy of the numerical integration and 
to test the correctness of a cosmological simulations. 

The evolution of this kind of structure in a 
static background, has been studied since the 

earlie st numerical N-body , 

1964t Ivan Albadal. Il982t lAarseth et al.l. llQSS : 



David fc Theunsl. Il989l: iTheuns fc Davidl. Il990t iBert 



200d: iHeggie k Hutl. '2003': Boilv fc Athanassoula'.~'200a 
Joyce. Marcos fc Svlos Lab ini. .2009a : Syl os Lab ini, 
2OI2I) . Briefly, driven by its own mean field, it collapses 
in a time scale « ^/Gp , forming a virialized quasi- 
stationary state (QSS), i.e. a collision-less configuration 
that is not in thermodynamical equilibrium but that 
has a finite lifetime fixed by the coUisional relaxation 
time-scale 



T2 ~ kNTc 



(1) 



where N is the system particles number and k is a nu- 
merical f actor that takes i nto a. c count the "Coulomb log- 
arithm" (ICh andrasekha?. ' 19431; 'Theis fc S purzeml . Il999t 
Diemand et al. . 2004: GabrieUi et al . 2010). 



Joyce fc Svlos Labinil (|2013[ ) have shown that an iso 



lated over-density in a cosmological simulation should 
in principle evolve, in physical coordinates, just as if it 
were in a non-expanding space with open boundary con- 
ditions. This corresponds to the fact that such a structure 
is in the stable clustering regime. Namely, it evolves as 
it were truly isolated from the rest of universe. While 
the stable clustering regime is assumed to approximately 
describe the e volution of a s tructure in the strongly non- 
linear regime ()Peeblesl . [l980l) . in the case of a single over- 
density with periodic boundary conditions and negligible 
finite size effects it must be precisely followed: the sta- 
ble clustering regime thus corresponds to a QSS for a 
structure in a static background. 

Indeed, in a cosmological simulation, where periodic 
boundary conditions are used, a single over-density can 
be considered isolated if finite size effects are negligible. 
Such mass distribution is isolated but for the interaction 
with its "copies" included in the infinite system over which 
the force is summed. It is easy to show that the copies 
contribute to the force with terms that represent pertur- 
bations suppressed by positive powers of Rq/L, where Rq 
is th e over-density s ize and L is the side of the p eriodic 
box (|Gabrielli et al.l . l200fil: Ijovce fc Svlos Labinil . l2013h . 
When the latter condition is satisfied, the only difference 
in the evolution between the expanding and the static case 
can arise from the smoothing of the gravitational force. 

Beyond the comparison of the expanding case with the 
template given by the static background evolution, one 



may consider the the analysis of the LI equation. In fact, 
Joyce fc Svlos Labinil (2013 ) found that is a useful diag- 
nostic to control the accuracy of cosmological simulations 
when the evolution of a single over-density is considered. 
In this case, the breakdown of the LI test occurs when 
the structure leaves the stable clustering regime: this is 
due to the difficulty of numerically integrating the system 
evolution with sufficient accuracy because hard two-body 
collisions take place. 



We consider in this paper the evolution of an ini- 
tially spherical and uniform cloud of particles with a 
velocity dispersion such that the initial virial ratio is 
in the range 60 G [~1jO]- When the system is warm 
enough, i.e. 5o ~ — 1, it approximately maintains its orig- 
inal size and shape during and after the collapse. Instead, 
when the system is cold enough, i.e. 60 ~ 0; it under- 
goes to a stronger contraction so that its size is consid- 
erably reduced during the collapse forming thereafter a 
virial ized structure with a density profile decaying as 
(IJovce. Marcos fc Sylos Labini l2009al : ISvlos Labini 
I2OI2 ). Our goal is to study the collapse of such a single 
cloud of particles in an expanding universe: by monitoring 
the LI equation and the stability of the clustering (in real 
space) we may understand to which accuracy the cosmic 
energy equation must be satisfied, in the strongly non lin- 
ear regime, so that the simulation still provides with phys- 
ically meaningful results. In particular we will study the 
accuracy of the numerical integration by changing the nu- 
merical and physical parameters that characterise a sim- 
ulation. 



The paper is organised as follows. In SectH] we re- 
call the main results of the evolution of an isolated over- 
density i n a static background . Then we summarise the 
results of lJovce fc Svlos Labinil (|2013t ) concerning the ap- 
proximate equivalence of an isolated cloud of particles in 
a static background with open boundary conditions and 
the same cloud when it is placed alone in a box of an in- 
finite periodic system with space expansion — i.e., when 
the simulation is run in cosmological conditions. In Sectl3] 
we discuss the preparation of the initial conditions, the 
choice of the parameters that control the numerical in- 
tegration and we briefiy recall the LI test. The compar- 
ison between the evolution of an isolated over-density in 
a static and expanding background is presented in Sec|H 
A series of systematic tests to determine the role of the 
integration parameters such as the gravitational force soft- 
ening length and the accuracy of the time-stepping is pre- 
sented in SectlSJ Then in Sectl6]we discuss several tests 
performed by changing the physical parameters of the ini- 
tial conditions, such as the number of particles, the sys- 
tem size and the velocity dispersion. Finally in Sect (7] we 
discuss the principal results and then we draw our main 
conclusions. 
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2. Isolated over-density in a cosmological and a 
static background 

Firstly we briefly recall the main features of the collapse 
of an isolated over-density in a static background with 
open boundary conditions (hereafter STAT case). Then 
we discuss the equivalence of the evolution of an isolated 
structure in STAT case and in an expanding background 
with periodic boundary conditions (EXP case). 

2.1. Isolated over-densities in a static background 

In the STAT case, the global collapse of an isolated and 
spherical cloud of self-gravitating particles, of density p, 
is determined by its mean gravitational field and it is 
characterised by the time scale 



37r 
32Gp 



(2) 



When the initial velocity dispersion is such that the virial 
ratio is 6o = '^Kq/Wq w — 1 (where Kq/Wq is the system 
initial kinetic/potential energy) the collapse consists in a 
few oscillations that drive the system to a virialized QSS: 
this maintains appro ximately the size Rp and shape of the 
initial configuration ( Svlos Labinil . 2012h . 

On the other hand, when the initial velocity disper- 
sion is small enough (or zero), the system undergoes to 
a stronger contraction reducing its size by a large fac- 
tor. In this case the features of the virialized QSS de- 
pend explicitly on the number of system particles (at 
constant density): for instance, when particles are ini- 
tially Poisson distributed, the intri nsic characteris t ic size 
in Eq |3l scales as » ([Aarseth et al. 



19881 



Jovce. Marcos fc Svlos Labinil . l2Q09aD . Other noticeable 
features are: (i) during the collapse a fraction of particles 
gain enough kinetic energy to escape from the system; (ii) 
the density profile after the collapse is characterised by 
the power law decay described by 



n(r) = 



1 + (r/re)4 ' 



(3) 



(that we denominate the quasi equilibrium profile — 
QEP) where Vc and ric are two constants. 

The key to understand the formation of this density 
profile is represented by the physics of the ejection mech- 
anism. The main characteristics of the virialized QSS that 
can be framed in a simple physical model introduced by 
Svlos Labinil (|2012h are 



n{r) ~ r , 
, , n(r) r^/^ 

$ r = -tA r 



r.-5/2 



(4) 
(5) 

(6) 



where Vr is the radial velocity and ^(r) is the so-called 
pseudo phase-space density. Note that, once the QSS is 
formed, its potential and kinetic energy reach a constant 



value: this fact represents one of the diagnostics to numeri- 
cally detect the quasi stationarity of t he particle configura- 
tion djovce. Marcos fc Svlos Labini . 2009a : Svlos Labini 
I2OI2I : I.Tovce fc Svlos Labinil . l2013h . 



2.2. Equivalence of the evolution in a static and in an 
expanding background 

Jovce &: Svlos Labinil ( 2013 ) have shown that when a cos- 
mological code is used to simulate an isolated structure (of 
size i?o) in an expanding universe, it should reproduce the 
same results, in physical coordinates, as that obtained for 
the same structure in open boundary conditions without 
expansion. Let us bricfiy recall in which limit this equiv- 
alence is verified. 

Dissipation-less cosmol ogical N-body s imulations solve 
numerically the equations (IPeeblesL [l980h 



2H- 



dt 



^ copie 

= rCm > 



copies j—l,..N 

E ■ 



(7) 



where x^ are the comoving particle coordinates of the 
i = 1, N particles of equal mass m, enclosed in a cubic 
box of side L, and subject to periodic boundary conditions, 
a{t) is the appropriate scale factor for the cosmology con- 
sidered, and H{t) = d/a is the Hubble constant. Note 
that the force on a particle is that due to the — 1 other 
particles and to all their copies. The leading non-zero 
term of the force due to the copies is a repulsive "dipole" ; 
by neglecting the quadrupole and higher multipole terms, 
that are convergent and suppressed by positive powers 
of (Rq/L) compared to the dipo le term ( Gabrielli et all . 
20061 : Jovce fc Svlos Labinil . 2013), we can rewrite EqlTjas 



dt^ 



2H- 



dt 



j=l,..N 
3^1 



(8) 



where the sum is now only over the A^ — 1 particles in 
the simulation box. Assuming for simplicity an Einstein- 
de Sitter (EdS) cosmology, for which 

d AnGpo 
a ' 

EqlTjmay be written, in physical coordinates = a{t)xi 
simply as 



(9) 



dt^ 



-Gm 



j=l...N 

E ■ 



(10) 



These are the equations of motion of A^ purely self- 
gravitating particles. Therefore, up to finite size correc- 
tions that vanish in the limit Rq/L — )■ 0, the equations of 
motion (in physical coordinates) in an expanding universe 
are the same of those describing the evolution of the same 
finite structure in a static background with open boundary 
conditions. 

We stress that EqlTUj represents an important approx- 
imation that justifies the analogy between the same cloud 
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of self-gravitating particles in STAT and EXP conditions. 
This is not a trivial approximation as it consists in replac- 
ing the sum in Eq|7]on a infinite periodic system, i.e. with 
infinite terms, with the sum on a finite system in Eo llOl 
The periodic sum can be written as a local sum on N par- 
ticles and a non-local term which takes into account the 
contribution of the infinite copies: this latter term can be 
neglected for the system considered, i.e. for Rq < L. Only 
in this situation the evolution of a single over-density in a 
periodic system, i.e., an infinite number of over-densities 
one in each of the infinite boxes of the periodic system, in 
an expanding background is equivalent to the evolution of 
a truly isolated over-density in a static background with 
open boundary conditions. 

Note that the simulations are run, as usual, in comov- 
ing coordinates. In what follows we present analyses both 
in comoving and in physical coordinates because the stable 
clustering regime should manifest itself in physical coordi- 
nates, i.e. the stru cture re mains almost invari ant only in 
these coordinates ( Jovce fc Svlos Labinil . 2013 ). 



3. The simulations 

3.1. Initial conditions 

We have considered a system with total mass M in a 
spherical volume V — Att/SRq: this can be discretized by 
using N randomly distributed particles of mass m such 
that 



M ZNm 



V 



— const. 



(11) 



We have used 10^ < TV < lO'^ (series S), TV = lO" (series 
M) and 10"* < TV < 10^ (series L). Details of the simula- 
tions are reported in Table [1] all length scales are given 
in units of the box side L. Note that the results we report 
require only choice of units for length and energy: for the 
former we will take units defined by L = 1, and for the 
latter units in which Wo = 3GMV(5i?o) = -1, I.e. m 
which the initial continuum limit gravitational potential 
energy is minus unity. 

Note that the system radius i?o is chosen to be i?o < 
L/4 so to minimise, as discussed above, the effect of the 
periodic boundary conditions. In particular, the distance 
between any pair of structure particle i,j is less than that 
separating i from any particle of the copies in the infinite 
periodic system. 

The spherical cloud corresponds, in a cosmological sim- 
ulation, to an over-density with initial mass density con- 
trast 



P_ 

Po 



3L3 
47ri?3 



(12) 



where p = 3M/{4ttRq) is the cloud density and po — 
M/L^ is the mean mass density of the universe (see 
Tabl21). Particles mass is such that po is equal to the crit- 
ical density. 



The scale factor in an EdS universe obeys to (jPeebled. 

19801) 



a{t) =ao[ — 



2/3 



(13) 



where to = 1/\/67tGpo. Assuining gp = 1 from Eqs lT3T^ 
we find ( Jovce fc Svlos Labinil . 12013 ) 



t-to = [a'/' - l]to = :^VS [0^/2 - i] 



(14) 



where Tc is given by Eql21 Note that from E cfHl we find 
that only for very high values of the over-density, i.e. S > 
10^, at moderate expansion factors a « 10, {t — tQ)/Tc is of 
the order of a thousand i.e. of the order of T2 for N = 10"^ 
(see EqlH). 

3.2. Numerical integration parameters 

We have used the publicly available tree-code GADGET 



(jSpringel et al.l . 120011 : ISpringell 120051 ). In this code one 



may choose the values of a number of parameters and 
we have followed the prescriptions of the user guide for 
GADGET-2 to fix the basic ones. 

The gravitational potential used by GADGET has the ex- 
act Newtonian potential above a separation of 2.8e while 
for smaller separations the potential first derivative van- 
ishes H. The different values of e used are reported in 
TabU] Note that we have chosen a softening length that 
is fixed in comoving coordinates and thus it varies in phys- 
ical ones. 

(l2003l) : 



Further we have set, following iPower et al 



ISpringel ( 2005 ) a Barnes- Hut opening criterion with an 
opening angle of the tree algorithm 6 = 0.7 for the first 
force computation and a dynamical updating criterion 
thereafter The accuracy of the relative cell-opening 



criterion (i.e., the parameter ErrTolForceAcc) □ is set 
to (j) — 5 X 10~^. Finally, we have chosen the time- 
stepping criterion of GADGET where the value of the 
parameter controlling its accuracy rj (i.e., the parameter 
ErrTolIntAccuracy) is reported in Tab|TJ 



3.3. Basic definitions 

Let us introduce some basic definitions that will be useful 
in what follows. The kinetic energy in physical coordinates 
is 

(15) 



and the potential energy in physical coordinates is 



^ see |http: // www. mpa-g arching.mpg ■ de/gadget /users-guide .pdf "I 

^ the value of e is the value of the parameter with this name 



in the code 
^ but 6* = 0.1 in the simulation M4b 
^ but (h = 5 X 10~* in the simulation M4c 
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where 



$(r,) = -G X! "^jffdJ'i - ""j 



(17) 



In order to test whether the LI equation is numerically 
satisfied during the time integration, one may measure the 
deviation of A{a) from 1. In particular, we will determine 
the expansion factor aLi such that 



and where q{r ) is the exact GADGET pair potential 
( Smingei [20051) . The virial ratio is defined as 



\A{aLi) - 1| < OT for a < gli , 



(27) 



5=^ 



(18) 



Note that the initial virial ratio is bo — 2Kq/Wq, where 
Kq , Wo are respectively the initial kinetic and potential 
energy. We define the potential energy in comoving coor- 
dinates 



Wc = ^ ^ mj$(xj) = Wp X a 

j 

and the GADGET gravitational potential energy 

copies 



- mj$(x, 



(19) 



(20) 



i.e., the two-body potential is calculated by using periodic 
boundary conditions. As mentioned above, in the limit 
Ro ^ L, we have Wg « Wc- 

The peculiar kinetic energy is 



where 

^ ' dt dt ^ ' 



(21) 



(22) 



and larger than 0.1 for a > aLi- This allows us to quan- 
titatively determine the range of redshifts in which the 
LI equation is well approximated in the simulations, i.e. 
in which the LI test is satisfied. Note that at very early 
times, when the initial velocities are close (or equal) to 
zero the denominator of Eal26l approaches zero and thus 
A{a) is difficult to be precisely determined. Otherwise for 
long enough times there is not any particular numerical 
problem in the determination of Ea l26l 

4. Comparison between the evolution in a static 
and expanding background 

We first compare the evolution of an isolated and uniform 
over-density in an expanding background with periodic 
boundary conditions (EXP) with the evolution of the same 
system in a static background and with open boundary 
conditions (STAT). We consider here the simulation M4a 
(see Tabll]), which we have chosen as a reference for the 
reasons illustrated below. Results for the other simulations 
are discussed in the next sections. 

Note that the difficulty of integrating the cold collapse 
occurs only around the strongest collapse phase of the 
system, i.e. at t w Tc, a time scale which is much smaller 
than any relevant value of the expansion factor considered 
in the what follows. 



If the system is viriahzed in a region of size Rq <^ L then 4.1. Short time evolution 



dvj 



GM 



dt J Ro 
In this approximation 
dvi 



Fig ID (upper panel) shows the potential and kinetic energy 

(23) in physical coordinates (normalised to the initial potential 
energy Wq). As for the STAT case the potential (kinetic) 
energy decreases (increases) to reach its minimum (maxi- 
mum) at i ss Tc, i.e. at the strongest phase of the collapse. 

(24) For t > Tc the bounded structure is stationary, i.e. Wp w 



so that the kinetic energy in physical coordinates (Eg fTSt 
is well approximated by the peculiar kinetic energy 
(EqlHI), i.e. Kp « Kc. 

3.4. The Layzer-lrvine equation 

Let us briefly recall the expression of the LI equation. By 
defining the total peculiar energy as E — {Kc + Wp) the 
LI equation can be written as 



-Kr 



d{aE) 
da 

Integrating Eql^Dwe find 

^^^^ ^ \aE{a) - aoE{ao)\ ^ ^ 
lao Kc{a)da 



(25) 



(26) 



const., and it satisfies the virial equilibrium (FigH] — bot- 
tom panel). A fraction of about /p « 25% of the particles 
gain kinetic energy so that, after the collapse, they are 
ejected from the system moving as free particles there- 
after. 

Energy conservation for the STAT case can easily be 
controlled to a precision better than « 0.01% (see, e.g., 
lAarsethl ( 2003f ) and references therein). From the LI test 
(FigH] — bottom panel) in the EXP case we find A{t) ~ 1 
with relatively large fluctuations up to the maximum con- 
traction phase. However, these fluctuations are damped 
after the collapse: we thus conclude that in this range of 
time the system is still in the regime a < gli (see Eq[77|). 

The density profile (Figl2] — upper panel) is well fit- 
ted by EqO both in the EXP and STAT case. The ra- 
dial velocity dispersion (Figl2] — bottom panel) shows 
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Fig. 1. Upper panel: Energy in physical coordinates {Kp 
Eg [island Wp Eg fTS)) normalised to the initial total energy 
as a function of time (computed by using Ea ll4l for the 
EXP case), together with the corresponding behaviour of 
the same system in the STAT case. Bottom panel: virial 
ratio \b(t)\, together with the LI test (Eql^H) for the EXP 
case. 



Fig. 2. Comparison of the density (upper panel) and ve- 
locity (bottom panel) profiles of the same system evolved 
in an expanding background (EXP at a = 1.36 — black 
circles) and in a static background with open boundary 
conditions (STAT at t = 5tc — red squares). The best fit 
with Eql3] (QEP) and a line with slope —1 are shown as 
references (see EqlS|). 



the Keplerian behaviour described by Eq|5] Therefore the 
pseudo phase-space density obeys to Eq|6l i.e. it decays as 

4.2. Long time evolution 

The situation discussed in Sect 14 . II drastically changes at 
longer times in the EXP simulation. Indeed, already for 
moderate values of the expansion factor, the evolution 
of the system rapidly deviates from the short-time be- 
haviour, and correspondingly the properties of the QSS 
become different from Eqs|4]|6l 

We find that Wp « const., i.e. the structure is in a QSS 
and in the stable clustering regime, only for a < a^i w 6 
(see FigH]). Instead, for a > qli we find that \Wp\ ^ 1/a 
(and thus Wc « const). Thus the breakdown of the LI 
test corresponds to a departure from the stable clustering 
regime. Indeed, the density profile in comoving coordi- 
nates for a > qli stabilises to an almost "asymptotic" 
shape which does not change anymore at small radii and 
slowly changes at large ones (see FigH]). This is charac- 
terised by an approximate decay extending over about 
three decades, i.e. from ^ IQ^^L to L. 

Note that the softening length is fixed in comoving co- 
ordinates and thus in physical coordinates it increases as 
Ep = as. When the structure is in the stable clustering 
regime, its size in physical coordinates Rp is constant. 



and thus Rp/Sp oc l/a: if the stable clustering regime 
persists long enough then the size of the structure will 
become of the order of, or even smaller than, the soften- 
ing length. However, we noticed that beyond a certain 
time the structure becomes stable in comoving coordi- 
nate and thus its size in physical coordinates scales as 
R* — aR*^, where i?*^ is the constant size in comoving co- 
ordinates. Therefore for a > a^i one approximately finds 
that Rp/Sp « Rcc/^ ~ const. 

The length scale R*^ can be estimated, for instance, by 
measuring the radius i?(0.4) (i?(0.8)) including the 40% 
(80%) of the mass of the virialized structure. As shown 
in inset panel of FigHJ for a > qli, i?(0.4) and i?(0.8) 
remain respectively larger than 10 and 100 times the scale 
beyond which the force is exactly Newtonian, i.e. 2.8£. In 
this condition the virial ratio fluctuates around -1 at all 
times (see Fig 13]). This implies that the Newtonian value 
of the mean field potential is well approximated by the 
smoothed potential and thus the softening length does not 
perturb the mean field dynamics. (Note that this result 
is confirmed by the behaviours measured by using other 
values of the softening length, as well as by an analysis of 
the gravitational potential energy in function of e — see 
discussion in Sect lS.ip . 

Coherently with the behaviour of the energy, the den- 
sity profile in physical coordinates (see Figl5]) does not 
change for a < a^j, and it is well fitted by EqISj this 
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Fig. 3. Upper panel: Absolute value of the potential en- 
ergy in physical coordinates Wp{a) (Eg fTOI) as a function 
of the scale factor. The approximate value of oz,/, sepa- 
rating the regime where Wp ~ const, from the one where 
\Wp\ ^ a^^, is reported as reference. Bottom panel: LI 
test (Eg 125)) and virial ratio (absolute value) for particles 
with negative energy. 



Fig. 4. Evolution of the density profile in comoving co- 
ordinates. Two reference lines with slope -4 and -3 are 
reported as a reference. The arrow shows the direction 
of the density profile evolution. Inset panel: behaviours 
in function of the scale factor of the length scale i?(0.4) 
(i?(0.8)), including the 40% (80%)of the mass of the viri- 
alized structure, normalised to 2.8e. 



situation corresponds to the QSS in the stable clustering 
regime. Note that, as time passes, the spatial extension 
of the halo, i.e. the range of radii where n{r) ^ r^^, 
increases. Indeed, as time passes particles with energy 
Bp < may go asyin ptotically far from the structure core 
(jSvlos Labiiiil . l2012l) . 

On the other hand, for a > gli, we find that n{r) 



continues to evolve: at large radii it decays as 



and 



at small radii the density does not flatten as for the STAT 
case. A iNavarro. Frenk fc White! (|l997l) (NFW) profile 



n{r) 



(r/r,)(l + (r/r,))2 



(28) 



gives a reasonable fit for a > uli (bottom panels of Fig[5]): 
in particular the large scale tail shows a well-defined 
decay. Therefore only when stability in physical coordi- 
nates is observed, the density profile is consistent with 
the expected one, i.e. EqlSj when stability in physical 
coordinates is broken then the density profile changes its 
functional behaviour approaching the NFW one. 

Similarly, the radial velocity dispersion profile in phys- 
ical coordinates (FigE]) shows an approximated Keplerian 
decay at large scale as in the STAT case. However, while 
for a < uli the different curves collapse on each other, 
implying that the virialized structure is in the stable clus- 
tering regime, for a > uli the velocity profile marks as 
well a breaking of the stable clustering regime. In par- 
ticular, the slope remains —1, i.e. halo particles are still 



orbiting in a central gravitational potential but with an 
amplitude that decreases with time. This is consistent with 
the fact that the potential energy in physical coordinates 
Wp increases with the scale factor for a > gli (as Wc « 
const.): as time passes R* ex. a and thus Wp increases pro- 
portionally to — 1/a. Thus for a > uli there is an energy 
injection into the system which causes an increase of the 
total energy: this occurs slowly enough that the system 
can rearrange its phase space configuration to satisfy the 
virial theorem (see the bottom panel of Fig [31). Therefore 
the energy increase corresponds to decrease of the kinetic 
energy: as 6 « — 1 we have AE = —AKp. Note the phys- 
ical size of the structure increases almost simultaneously 
with the violation of the LI equation: this implies that the 
energy injection is spurious as otherwise the LI test should 
be verified. 

The decrease of the kinetic energy corresponds to the 
fact that high velocity particles are slowed down. This is 
shown in FiglT] (upper panel) where it is reported the evo- 
lution of the probability distribution function (PDF) f{v) 
of the absolute value of the velocity. While for a < a^i we 
find that f{v) remains almost stable, for a > a^i high ve- 
locity particles are slowed down and correspondingly the 
high velocity tail of f{v) is significantly reduced. However 
the shape of f{v) remains almost the same for a > oli- 
The failure of the LI test and consequent breakdown of the 
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Fig. 5. Density profile in physical coordinates for the sim- 
ulation M4a at different scale factors a (see labels) to- 
gether with the density profile for the STAT case. The 
arrow shows the direction of the density profile evolution 
in the cosmological case. Further it is shown the best fit 
given by Eql3] (QEP) and a reference line with slope -3. 




r 



Fig. 6. Radial velocity dispersion profile in physical co- 
ordinates at different scale factors a (see labels) together 
with the profile for the STAT case. The arrow shows the 
direction of the density profile evolution in the cosmolog- 
ical case. A line with slope r^^ (see EqO is reported as 
reference. 




log{v) 

Fig. 7. Evolution of the probability distribution function 
(PDF) of the absolute value of the velocity at different 
scale factors (see labels). The arrow shows the direction 
of the evolution. 

stable clustering regime is thus induced by the difficulty 
of integrating high velocity particles. 

High velocities are typically associated with close scat- 
terings and the limited ability to resolve hard scattering 
between close neighbours, through which particles can in- 
crease their velocity, can also be traced by looking at the 
evolution of the nearest neighbour distribution uj{x) in co- 
moving coordinates (see FigjT]): indeed, one may note that 
while for a < a^j lo{x) shows a peak that is shifted at 
smaller and smaller scale, corresponding to the fact that 
the structure is stable in physical coordinates, for a > Uli 
uj{x) does not change anymore reaching an "asymptotic" 
shape. This is a clear effect of the limited small scale res- 
olution of the numerical code. 

Finally we note that from the behaviours of the density 
profile and of the radial velocity dispersion, we can easily 
find that the pseudo phase-space density, for a < uli and 
for r > Tc, scales as $(r) ^ r~^/^, i.e., the behaviour 
expected for the QEP (see EqlH]). On the other hand, for 
a > aLi: at large radii, we find $(r) ~ x r^/^ ~ r~^-^. 

4.3. Summary 

Let us summarise the situation in the two regimes that we 
have identified., i.e. a < a^i and a > a^/. 

1. For a < aLi the structure is in a QSS (i.e., \Wp\ « 
const.) and the density and velocity profiles are de- 
scribed by the collision-less equilibrium behaviours 
(EqsHlini). This situation corresponds to the stable 
clustering regime: the physical size and shape of the 
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Note that a^j k, 6 corresponds to a time (see Eq lTil) 
(tbreak — ta)/Tc ~ 50 that is much smaller than the col- 
lisional relaxation time scale (see EqU]) T2 ~ W*Tc for 
a system with N = 10'* particles. This implies that the 
break-down of the QSS is not related to the expected evap- 
oration due to two-body collisions for t > t^. 

Finally we note that the breakdown of the QSS, if 
physically originated, should not be associated with the 
breakdown of the LI test as we find. This shows that such 
a breakdown is induced by spurious reasons that we iden- 
tify in the role of the inaccurate integration of high veloc- 
ities particles and thus of hard scatterings. We will return 
on this important point in Sect l5.1| 



10 

x/{2M) 

Fig. 8. Evolution of the nearest neighbour distribution in 
comoving coordinates at different scale factors (see labels) . 
The arrow shows the direction of the evolution. The initial 
distribution (IC) is also shown for comparison. 

structure remains invariant during time evolution, 
modulo a marginal spatial extension of the power-law 
tail due to high velocity, but bounded, particles. 
2. For a > a^j the structure is (almost) frozen in co- 
moving coordinates and therefore its size in physical 
coordinates increases and its potential energy tends 
to zero as Wp oc —1/a. The breakdown of the IL 
test marks a problem in the numerical integration of 
the particles equations of motion in this regime: this 
occurs because there is an injection of energy which 
causes the expansion of the structure in physical space. 
Correspondingly the density profile develops a 
tail, approaching a NFW profile. The velocity dis- 
persion profile still decays in a Keplerian way but with 
an amplitude that decreases with time equilibrating 
the smaller potential energy (in absolute value) of the 
structure, so that the virial ratio remains & ss — 1. 

The decrease of the system kinetic energy for a > gli 
corresponds to the fact that high velocity particles are 
slowed down and the PDF of the particles velocity f{v) 
shifts toward smaller velocity values. This is indeed the 
consequence of the inaccurate integration of two-body col- 
lisions: it becomes more difficult to numerically resolve 
high velocity particles. As shown by the behaviour of the 
nearest neighbour distribution uj{x) these high velocity 
particles are also close neighbours. Indeed, ui{x) does not 
evolve anymore for a > ol/, as it should if stable cluster- 
ing were satisfied. We thus conclude that the inaccurate 
integration of hard collisions induces the change of regime 
at uli- 



5. Role of numerical integration parameters 

In this section we study the effect of varying some of the 
most relevant numerical parameters that characterise the 
numerical integration: the softening length e, the param- 
eter that controls the time-stepping accuracy 77, the open- 
ing angle 9 of the tree algorithm and the accuracy of the 
relative cell-opening criterion 00. 



5.1. Softening length 

The softening length |^ has been varied in a range such 
that: (i) e < Rq, i.e. it is smaller than the initial size of 
the structure so that the gravitational mean field force 
is well approximated at least at early times and (ii) the 
lower limit to e is due to numerical constraints: e > 10~^ 
otherwise the integration time rapidly becomes too long 
and the simulation becomes unfeasible. 

The potential energy in physical coordinates Wp and 
the IL test (FiglH]) for simulations of the same system with 
N = 10* particles and with different values of e (see Tab[T]) 
— changed by a factor 10* — show behaviours coherent 
with those discussed in Sect 14. 21 Namely the two regimes 
found above, corresponding to Wp « const, for a < uli 



and \Wp\ ^ a^^ for a > uli, still characterise the temporal 
evolution. However, the transition between them depends 
on the softening length, i.e., a^i ~ aLi{s). 

In particular, we find that a^i w 6 for e € [3.7 x 
10~^,3.7 X 10"*], while out side this range a_Lj ~ 1 . This 
result is in agreement with Joyce fc Svlos Labinil ( 2013 ) 
who concluded that at most in that narrow range of e 
can the behaviour required, a QSS in the stable clustering 
regime, be reproduced by the N-body method. 

In order to understand the effect of a softening length 
which is too large with respect to the size of the structure, 
we have computed the gravitational potential energy of 
the structure with density profile described by Eq|31 in 
the case in which the pair potential is the GADGET one with 
different values e. We find that the gravitational potential 
energy is very well approximated by the Newtonian value 



^ i.e., the parameter denominated ErrTolForceAcc in the 
code 

® recall that e is given in units of the box side L 
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Fig. 9. Behaviour of A{a) (Eql^H) for simulations with 
different e (see TabUfor details). 



when i?(0.4)/(2£) « 10 and/or i?(0.8)/(2e) « 100, where 
i?(0.4) and i?(0.8) are the radii in comoving coordinates 
including respectively the 40% and 80% of the mass of the 
virialized structure. 

The behaviours of i?(0.8)/(2e) as a function of the ex- 
pansion factor and for different values of the softening 
length are reported in FiglTU](see the bottom panel). For 
e > 3.7 X IQ-^ we have that i?(0.4)/(2.8£) < 10: in this 
situation the effective mean field force due to all particles 
is different from Newtonian, i.e. the smoothening length is 
too large with respect to the size of the system. The effect 
of a too large softening thus corresponds to a deviation of 
the virial ratio from -1 (see the upper panel of Fig fTU)) . 

Instead, for e < 3.7 x 10^^ we find that, be- 
cause the system remains stable in comoving coordinates 
i?(0.8)/(2.8e) is larger than 100 for all values of the 
scale factor considered, and thus the virial ratio fluctu- 
ates around -1 at all times. Therefore we conclude that 
the relevant characteristic size of the structure remains 
large enough with respect to the force softening and that 
in this situation the Newtonian value of the mean field 
potential is well approximated by the smoothed potential. 

The behaviour of the density profile is also coherent 
with what we have already found in Sect|4l for a < a^i, 
the density profile obeys to Eql21 

Then, for a > a^j, we observe a breaking of the stable 
clustering regime, and the density profile rapidly changes 
shape. In particular, it firstly deviates at small radii, dis- 
playing a power-law behaviour, and then at large ones 
where the slope tends to —3. The whole density profile 
can be fitted by a NFW behaviour (see Fig fTTj) . 



cosmological N-body simulations 




Fig. 10. Upper panel: virial ratio as a function of the 
scale factor in simulations with different e. Bottom panel: 
Length scale i?(0.8), including the 80% of the mass of the 
virialized structure, normalised to 2.8e. 




r 



Fig. 11. Density profile in physical coordinates at a = 4.6 
for the simulations M4a and Mia. The best fit with the 
NFW and QEP profiles are also reported. 

The evolution of the the velocity PDF f{v) and the 
nearest neighbours distribution uj{x) for different values of 
e are coherent with the behaviours previously discussed. 
Namely, for large enough scale factors, i.e. a > 10, the sim- 
ulation which behaves better according to the LI test (i.e., 
£ = 3.7 X 10~^) shows a more extended tail of f{v) at high 
velocities than the others. This means a better ability of 
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Fig. 12. Behaviour of A{a) (Eg 125)) for simulations with 
different 77 and e (see labels) 



the code to follow high velocity particles. Correspondingly 
also llj{x) is peaked at smaller separations i.e. close parti- 
cles can be better resolved for that value of the softening 
length. Both these facts point toward a better spatial (i.e., 
close neighbours) and temporal (i.e., high velocity parti- 
cles) resolution when e = 3.7xl0~^. 

Similarly to the behaviour of the density profile, whose 
small scale properties depend on e the nearest neighbour 
distribution and its peak shown an e dependence implying 
that the small scale properties of the system are deter- 
mined by the spatial resolution used in the simulation. 

5.2. Time-stepping accuracy 

The accuracy of the time-stepping in the GADGET code is 
controlled by the parameter rj: in particular, the time- 
stepping is fixed by 

At^JJl, (29) 

where Aacc is the particle acceleration. It thus depends 
both on the softening length e and on the parameter rj, 
i.e. the spatial and temporal resolution are strongly entan- 
gled. In order to determine the role of rj we consider a set 
of simulations in which all physical {N, Rq, bo) and numer- 
ical parameters (e, 9, (p) are fixed and the time-stepping 
accuracy 77 is varied by a factor 10^ (see Tab[T]for details). 

From the LI test (see FiglT^ we find aLi = aLiiv)- 
Note that, as previously found, the LI test is verified in 
about the same range of scale factors where the poten- 
tial energy is roughly constant and the structure is in the 
stable clustering regime. 

The behaviour of a^i for different values of e and 77 is 
reported in Fig ll3l the larger is 77 and the sooner occurs 
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Fig. 13. Behaviour of a^/ as a function of 77 for different 
values of the softening length e. 
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Fig. 14. Density profile in physical coordinates at a = 2.5 
respectively for 77 = 0.0025 and 77 — 0.25. In both the cases 
the softening length is e = 3.7 x 10^^. 



(i) the breakdown of the LI test, (ii) the departure from 
the stable clustering regime, (iii) the deviation from the 
QSS with the profile given by EqlH] That is, at fixed e, the 
larger is 77 the smaller is a^/. In agreement with the results 
discussed above, we find that the NFW profile provides a 
very good fit of the density profile for a > a^j (see Fig fH)) . 

The behaviours of the velocity f{v) and of the nearest 
neighbours uj{x) distributions show that the smaller is 77 
and the better is the resolution of close neighbours and of 
high velocity particles. In addition, the high velocity tail 
of f{v) is substantially reduced in the regime a > uli, 
in agreement with the results obtained by changing the 
softening length. 



12 



Sylos Labini: Testing the accuracy of cosmological N-body simulations 



5.3. Other numerical parameters 

Results obtained by varying the opening angle 6 of the tree 
algorithm from 0.7 to 0.1 and the accuracy of the relative 
cell-opening criterion 0, i.e. the parameter denominated 
ErrTolForceAcc, from 5 x 10""^ to 5 x 10^^ do not show 
sensible differences. This fact gives us reasonable confi- 
dence that the results we have presented are quite well 
converged with respect to these parameters. 



5.4. Summary 

The two main parameters that control the spatial and 
temporal resolution of the numerical integration, namely 
the softening length e and the time-stepping accuracy tj, 
determine the range of time in which a simulation is well 
performing, i.e., the structure is a QSS and Wp ~ const., 
the evolution is in agreement with stable clustering and 
the LI test is reasonably well satisfied. Thus the scale 
factor at which the breakdown of the LI test occurs is 
such that gli — aLi{e,i^) . We note that most of the 
discussion in the literature has focused the attention on 



allowing a better sampling of the change of velocity during 
hard scatterings. 

Instead, the change of particle velocity is Aw At x 
Aacc so that for hard scatterings Aacc cx we have Aw oc 
yjrjje: the over-estimation of the velocity change increases 
when £ gets smaller. Therefore we may conclude that the 
measured dependence of a^i on e, 77 and N is coherent 
with the interpretation that there is a numerical problem 
in the integration of hard scatterings. CoUisionality for the 
smaller softenings we use manifests itself most clearly in 
the breakdown of the LI test that represents a useful tool 
to study these effects. 

It is interesting to note that collisional effects in cos- 
mological simulations have been documented in previous 
studies by co nsidering the c ompa rison of different codes. 



In particular, iKnebe et al 



the ro le of the softening length ( see e.g. ISplinter et al. 



(Il988l): iHeitmann et alT (I2OO5I) : iRomeo et al.l (l2008l): 



dlooa) very explicitly demon- 
strates the effects of poor integration of hard two body 
collisions at small softenings. They showed that two-body 
scattering have effects which can be greatly amplified by 
inaccuracies of time integration and that the effects of 
such poorly integrated hard collisions is to lead to an ar- 
tificial injection of energy into the system which causes 
an increasing its size. These results are indeed fully com- 
patible with our results and interpretati on (see also the 



Jovce. Marcos fc Baertschigeij (|2009t ): IPower et all discussion in lJovce fc Svlos Labin ? (|2013l )). 



( 2003[) ): here we have found that a crucial parameter is 
also represented by the time-stepping accuracy. 

Coherently with these results we found that the ve- 
locity PDF f{v) and nearest neighbour distribution ijj{x) 
depends on both the two parameters e and 77. In particular 
the high velocity tail of /(u) remains stable for a < a^i 
while there is a sensible deficit of high velocity particles 
for larger scale factors, due to the limited temporal resolu- 
tion. Similarly w(x) does not change anymore for a > a^/, 
implying the departure from the stable clustering regime. 
In this situation the peak of w(x) is determined by the 
values of both parameters £,rj. 

In summary, tests performed by varying the softening 
parameter show that at fixed £ < € the smaller is ry and 
the larger is a^/. Our interpretation of this behaviour is 
that at a given level of numerical accuracy, determined by 
the choice of e, ry, hard scattering between particles can 
be resolved only up to a certain scale factor. Finally we 
note that the results of the evolution of the density profile 
are coherent with the ones found in SectUl namely for 
a < o-Li this is well approximated by the QEP shape, 
while for a > a^i the profile shows the NFW behaviour. 

We have discussed that there is an artificial injection 
of energy into the system during time evolution. We now 
conjecture that this is due to an inaccurate integration of 
two-body scatterings. When an hard scattering occurs the 
acceleration can be huge for a very short time interval St: 
however if the time stepping is larger than 6t then the 
particle velocity will be overestimated by a large factor: 
this is a known numerica l problem discussed in detail by, 
e.g., Knebe et al. (j2000[ ). Indeed, when 77 decreases the 



amplitude of the time stepping (see Eal29)) gets shorter 



6. Role of the physical parameters 

In this section we first study the evolution of systems with 
same physical and numerical parameters but different par- 
ticles number. That is, we discretize the system with a dif- 
ferent number TV of particles taking fixed the total mass 
M = m X N , with m oc 1/N. Then we consider systems 
of different size Rq but with the number of particles, their 
mass and the numerical parameters [£,77,^,0] fixed: this 
corresponds to change only the value of S in EqfT^ Finally 
we change only the initial velocity dispersion and thus the 
initial virial ratio bg. 



6.1. Scaling with particles number 

When the number of particles is small enough, i.e. N < 
N* « lO'' we observe that during the collapse phase the 
LI test always fails as A{a) presents relatively large devi- 
ations from unity since a > 1 which persist at all times 
(see Fig fTS)) . This can be interpreted as due to the fact 
that a few hard collisions occurring in the collapse phase 
substantially perturb the integration of the equation of 
motions. Note that according to the results obtained in 
SectlS]the value of N* should in general depend on £,77. 

In simulations with N > N* « 10"* we find instead 
that A{a) w 1 soon after the collapse phase, with smaller 
fluctuations for larger N values. The scale factor signing 
the break of the LI test weakly depends on TV, i.e. uli = 
o-LiiN): indeed it is found to be almost independent on 
N. This latter result seems surprising as from EqU] the 
collisional relaxation time should increase with N so that 
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Fig. 15. Behaviour of A{a) (Eg 125)) for simulations with 
different N (see TabU]) Insert panel: Behaviour of the 
virial ratio for the same simulations 



one would naively expect that larger is N and larger is 
aLi- 

To interpret this behaviour we recall that the density 
of the core after the collapse, ric in Eq[31 scales with the 
number of particles as ric (x TV^: at fixed mass density, 
by increasing the number of particles, the density of the 
core increases as pc = ric x m (x N, given that particles 
mass s cales as m cx N^^ (jjovce. Marcos fc Svlos Labini . 
2009al) . From Eq[T]we thus find that T2 = kNt^ oc kVN, 
as Tc « 1/y/Gpc- Given that k is weakly dependent on 
TV, i.e. K cx l/log(iV), we thus find that T2 is almost in- 



dependent on A^, at least for the limited range of N we 
considered, i.e. N — 10^ 10^. 

In Fig ll6l it is shown the density profile in physical co- 
ordinates for simulations with N = 10^, 10^ particles: note 
that n(r) has been normalised by rec alling the results of 
Jovce. Marcos fc Svlos Labinil ( 2009a ) where it was shown 
that the two parameters in Eql3] scale as rc ~ N~^^-^ and 
Uc ~ iV^. As we found in the previous sections, the break- 
down of the stable clustering regime corresponds to a de- 
formation of the profile. In the case shown in FigfTBl while 
for TV = 10^ the simulation is still in the QSS regime, for 
the = 10'^ case the density profile is already well fitted 
by the NFW shape (see FigfTC)). 



6.2. Finite size effects 

By considering the behaviours of the gravitational poten- 
tial energy and of the LI test for simulations with different 
initial size Rq we can conclude that a^/ = aLj{Ro). In par- 
ticular, the smaller is Rq and the smaller is a^/. This fact 



Fig. 16. Density profile in physical coordinates for simula- 
tions with 10'^ and 10^ particles, at a = 4.9, renormalised 
as explained in the text at different values of the expansion 
factor. The best fit with EqUHl (NFW) and EqE] (QEP), 
respectively for N = lO'^ and N = 10^ are reported as 
reference. 



can be simply interpreted by considering that when the 
system size decreases, the amplitude of the over-density 
5 = 5{Rq) increases and so the value of the character- 
istic time scale Tc = Tc{Ro) (see Tabl2]). Thus, to com- 
pare the behaviours for different Rq we plot (see FigfTTl 
— bottom panels) Wp{a) and A{a) as a function of time 
in units, for each value of Rq, of the corresponding value 
of Tc ~ Tc{R()). The residual difference in the observed be- 
haviours is due to the effect of the different ratio Ra/L: 
when Rq > 1/lOL finite size corrections, i.e. the terms 
neglecte d in approx imating Eq|7] with Eq llOl become im- 
portant ( Jovce fc Svlos Labinil . I2OI3I) . 

The evolution of the density profile in physical coordi- 
nates for the systems with different Rq (see Fig lTSl) shows 
that they firstly reach the QEP behaviour; then the den- 



sity profile approaches the behaviour, at different val- 
ues of the scale factor for different Rq and in a coherent 
way with the behaviour of A{a). 

Finally we note that for long enough values of the scale 
factor the density profile for different Rq converges pre- 
cisely to the same behaviour characterised by a decay. 
As discussed in SectlSJ at fixed time-stepping accuracy, the 
small scale properties of the distributions are determined 
by the value of the softening length. Indeed, independently 
on the initial value of Rq, for long enough times, all curves 
collapse on each other. 
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Fig. 17. Behaviour of A{a) (Eg 125)) for simulations with 
different initial size Rq as a function of time (normalised 
to Tc{Ro)). 

6.3. Role of the initial velocity dispersion 

The initial conditions we have considered can be thought 
to represent an halo which virializes in a relatively short 




Fig. 18. Density profile in physical coordinates at a = 2.5 
for different values of Rq (see labels) . The vertical dashed 
line correspond to a x e. A reference line with slope —4 is 
plotted as reference. 



Fig. 19. Density profile in comoving coordinates at a = 70 
for the simulations with e = 3.7 x 10"'' and a different 
value of the initial radius Rq (see labels) . A reference line 
with slope —4 is plotted as reference. 



time and which then evolve in isolation from the rest of 
the mass of the universe. Indeed, when bo ^ —1 the halo 
is already almost virialized while when bo —?' the halo 
is far from such a configuration. However for bo £ [^1,0] 
the structure undergoes to a collapse, which can be less or 
more violent, that drives the system, in about the same 
time-scale Tc (see EqlJ), into a virialized quasi-stationary 
state that can be w ell described as a co llision-less equilib- 
rium configuration (jSvlos Labinil . |2012[ ) . 

FigEOl shows the behaviour of A{a) for simulations 
with different initial virial ratio bo G [—1,0], i.e. with 
different initial velocity dispersion. We find that uli — 
a-Liibo)- in particular, the smaller is bo the larger is uli- 
This behaviour is easily explained by considering that the 
closer the system is initially to the virial configuration and 
the less violent is t he dynamics that dr ives the system into 
the virialized QSS ( Svlos Labinil . |2012|) : the larger is bo at 
the higher is the density of the core reached after the col- 
lapse. 

We recall that, in the STAT case, the formation of 
the r~'^ is observed only when the collapse is violent 
enough that a fraction of the initial mass and energy is 
ej ected from the syst em at the formation of the QSS: 
in ISvlos Labini ( 2012 ') it was found that this occurs for 
&o > —1/2. In this situation the system evolution is very 
similar to the &o = case: at the beginning the density 
profile is well approximated by the QEP and then, for 
a > qli it develops a tail approaching the NEW shape 
also at small radii. Thus, although the 6o = is peculiar 
because the collapse is instantaneous in the continuous 
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Fig. 20. Behaviour of A{a) fEa l^^ for for simulations 
with different initial size bo. 

limit, it does not represent a pathological case and the 
tail is formed as long as the collapse is violent enough that 
part of the mass and energy is ejected from the virialized 
structure. 

For values of the initial virial ratio bo < —1/2 the 
virialized structure after the collapse does not present a 
density profile with a tail. However for long enough 
times the density profile approaches the 1/r^ behaviour for 
any value of bo. This shows that the 1 /r^ is formed through 
the dynamical process and that it is thus independent on 
the initial conditions. 

In Figl2T] it is shown the comparison of the density 
profile, in comoving coordinates, for simulations with dif- 
ferent initial virial ratio and at different expansion factors: 
this confirms that the smaller is bo and the slower is the 
evolution. In particular, the smaller bo the longer it takes 
to form the tail, a behaviour that is approached at 
long enough times for any value of bo- 

6.4. Summary 

In summary we found that qli = aLi{Ro, N,bo), i.e. the 
regime in which the system evolution is correctly numeri- 
cally integrated depends on the physical parameters of the 
initial conditions. 

For what concerns the dependence of a^j on N, we 
have noticed that: (i) for relatively small particles num- 
ber, i.e. N < N* « 10^, the system never reaches the 
quasi-stationary state in the stable clustering regime. We 
have interpreted this behaviour as due the fact that a few 
hard scatterings are enough to perturb the macroscopic 
dynamics, (ii) On the other hand for large enough parti- 
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Fig. 21. Comparison of the density profile in comoving 
coordinates for simulation with different initial virial ratio 
a, a ~ 70. 

cles number, i.e., N e [10"', 10^], we find that the deviation 
from the QSS slowly depends on N. This occurs because 
when initial conditions are cold, the system contracts more 
for larger N values, so that the virialized state has a core 
mass density pc oc N and thus the typical time scale for 
crossing is Tc cx 1/ \/N : in this situation the collisional 
time scale T2 (see EqlT]) is very weakly dependent on N. 

By decreasing the system size Ro we find that oli 
decreases as well: this is coherent with the fact that the 
smaller is Rq and the smaller is the average distance be- 
tween nearest neighbour particles in the core, and thus the 
higher the probability of hard collisions. However when we 
renormalize the behaviours to the proper characteristic 
time Tc ~ Tc{Ro) we find that the differences in the de- 
parture from the stable clustering regime become smaller: 
these residual differences are due to finite size effects. 

In addition we have found that for long enough times 
all systems with initially different _Ro converge to the same 
density profile characterised by a r~^ decay. This con- 
vergence is obtained because the numerical integration is 
badly performed, as shown by the LI test. 

Finally, by increasing the initial velocity dispersion, 
the breakdown scale factor qli is found to increase. This 
behaviour can be simply explained by considering that 
when 6o — ^ the collapse is much more violent than for 
bo — > —1; thus in the former case the core may reach 
higher densities sooner, then having a similar evolution to 
the cold case. In addition, when the initial velocity disper- 
sion is high enough, the system remains in the stable clus- 
tering regime longer and until the softening length, that 
increases in physical coordinates, becomes large enough to 
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perturb the evolution. At later times the system shows an 
evolution that is again similar to the cold case. 

7. Discussion and Conclusion 

In this paper we have tested a cosmological code (GADGET) 
by choosing very simple and idealised initial conditions. 
While we do not claim that this choice can be justified 
by any physical argument, i.e. the toy model does not 
represent any realistic physical system, we stress that in 
this way we have been able to make a series of controlled 
tests on the accuracy of a cosmological code. Given what 
we have learnt from this case, we discuss here some gen- 
eral ideas for the testing of a cosmological code in the 
case of more realistic initial conditions, as those predicted 
by standard models of galaxy formation. In order to con- 
sider separately these two different but, as we argue below, 
probably related issues, we first discuss the results of the 
simulations we have presented in this paper and then we 
briefiy comment on the more complicated problem that 
requires more tests and work to be properly explored. 



the dynamical point of v iew: indeed, as firstly noticed by 



Stiavelli fc BertinI (|1987l ). a density profile with a r tail 



7.1. The con 

code: a toy model 
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The evolution of an isolated over-density in a cosmological 
N-body simulation is approximately equivalent, in phys- 
ical coordinates, to the evolution o f the same structure in 
an op en system without expansion ( Joyce fc Svlos Labinil . 
2013r) . This equivalence is strictly verified in the limit of 



vanishing force smoothing and when the size of the over- 
density is much smaller than the size of the simulating 
box. When these conditions are satisfied, the evolution 
of the open system can be used as a template to com- 
pare the evolution in an expanding background with the 
aim of determining the accuracy of a cosmological sim- 
ulation. The open system evolution is such that it re- 
laxes to a virializ ed quasi sta.tionary state (QSS ) in a 
time Tc (see Eal^ (Ivan Albadal . ll982t Aarseth et all . ll986 ; 
Jovce. Marcos fc Svlos Labinil . l2009at ISvlos Labini 12012 ) 
and it evaporates, due to coUisional effects, in a time 
^2 ^ Tc (see Eq|T]). The QSS corresponds, in a cosmo- 
logical N-body simulation, to an isolated structure in the 
sta ble clustering regime. 



Jovce k, Svlos Labinil (|2Q13l ) considered the evolution 



of an isolated over-density initially close to virial equi- 
librium: for this reason during the collapse its size and 
density vary a little. In this paper we considered a more 
general case, namely an isolated over-density with initial 
virial ratio in the range 69— € [—1,0]. When 60 ^ ~1 
the initial conditions are warm enough that the collapse 
consists in a relatively small modification of system size. 
Instead when initial velocities are cold, i.e. 60 — ^ 0, the 
system largely contracts during the collapse, reducing con- 
siderably its size and shape and thus rapidly reaching, i n 
its central core, a very high density (|Svlos Labini |2012| ). 

Note that an initially cold, i.e. 60 = 0, and spherical 
over-density does not represent a pathological case from 



is a generic property of systems forme d by a violent relax- 
ation process. For instance, recently ISvlos LabiJ (|2012l 
I2OI3I) found that this profile is formed also in the col- 
lapse of an initially spherical and isolated cloud of particles 
with a small enough velocity dispersion, with a uniform 
or power-law density profiles. 

We found that, for 60 0, once the system has reached 
a virialized configuration, its evolution can be charac- 
terised by two different regimes: 

— (i) for scale factors a < qli the potential energy in 
physical coordinates Wp is approximately constant, 
the Layzer Irvine (LI) equation, describing the energy 
change in an expanding background, is well satisfied, 
the density profile is characterised by a r^^ tail, char- 
acteristic of the quasi equilibrium profile (QEP) in the 
open case (see Eq|3l), and the system is in the stable 
clustering regime, which means that it remains stable 
in physical coordinates. 

— (ii) On the other hand for a > gli we find that 
\Wp\ ^ a^^, there is a breakdown of the LI equation, 
the density profile becomes stable in comoving coor- 
dinates, which corresponds to the breakdown of the 
stable clustering regime, it develops a r^^ tail and it 
can reasonably be fitted by the NFW profile. 

In principle, the QSS has a lifetime T2 that diverges 
with the number N of system particles because of coUi- 
sional effects. However we found that the time scale cor- 
responding to a^i is (i) t^i < T2 and (ii) it depends on 
the numerical parameters of the simulation. Therefore this 
time scale is a numerical artifact: the relevant question is 
what determines it. 

In order to answer to this question, we have performed 
various series of simulations by varying the integration pa- 
rameters, such as the softening length of the gravitational 
force £ (which, in the simulations we have performed, is 
fixed in comoving coordinates during the evolution) and 
the accuracy of the time-stepping 77. We found that 

a-Li = aLi(e,rj). 

In particular, we found that there is a narrow range of e 
(at fixed rj) for which the simulation behaves "well" (i.e. 
regime (i) above) . When e is much smaller than the initial 
inter-particle distance l the integration is more difficult 
and aLi decreases: this fact can be interpreted as due the 
effect of hard collisions which require a high numerical ac- 
curacy to be correctly integrated. On the other hand, when 
£ is of the order of the characteristic scale of the virialized 
system the mean field potential is no longer Newtonian so 
that the system deviates from the virialized configuration 
and from the quasi equilibrium profile. 

That there is a problem in the accuracy of the nu- 
merical integration is confirmed by varying 77 at fixed e 
(with e < ^): in particular when rj decreases a^/ increases. 
This fact can be again interpreted in terms of a better 
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accuracy of the numerical integration and in particular a 
more accurate integration of high velocity particles, which 
require a smaller time-stepping to be properly followed. 
Both tests by varying e and 77 are thus consistent with 
the inter pretation outlined ab o ve, th at confirms the re- 
sults by IJoYce k. Svlos Labini ( 2013 ): two-body scatter- 
ings perturb the accuracy of the simulation evolution. 
Thus the temporal and spatial resolution of the simula- 
tion are entangled and both must be considered for the 
determination of the optimal values of the numerical pa- 
rameters. 

Note that codes with a dapti ve spatial resolu- 



tion (see e.g., Knebe et_alj (|2000l) : iTevssieil (|200l[ ): 



Miniati &: Colella (,2007. )') may perform much better than 
the fixed resolution case which we have considered in this 
paper. Indeed, it is clear such codes will, by construction, 
avoid respect. The tests we have discussed in this paper 
can be easily applied to other codes and a detailed study 
of codes with an adaptive mesh refinement is postponed 
to a forthcoming work. 

Because of the violation of the LI, for a > qli, the 
system expands in physical coordinates, thus increasing 
its potential energy: this is due to an injection of en- 
ergy into the system. From the tests performed by vary- 
ing the softening length and the time-stepping parameter 
we conclude it is the effects of poor integrated hard two 
body collisions that leads to an artificial injection of en- 
ergy into the system . This same conclusion was drawn by 
Knebe et al. ( 2000l) by considering a comparison of differ- 
ent N-body codes. In s umma ry we find, in agreement with 
Joyce fc Svlos Labinil ( 2013 ). that given a certain choice 
of numerical parameters, one may follow structure forma- 
tion in the non linear regime, for a limited time which 
depends on several physical properties of the system. The 
results obtained here clearly show the role of the numeri- 
cal (in) accuracy in the determination of the shape of the 
density profile structure. 

In order to further test the result that the numeri- 
cal integration becomes inaccurate because of coUisional 
effects, we have, at fixed numerical parameters, changed 
the physical parameters of the initial conditions. Namely 
we have varied the number of particles N (at fixed den- 
sity), the size of the structure Rq and the initial velocity 
dispersion (and thus the initial virial ratio bo). As result 
we found that 

aLi = aL/(iV, i?o,6o) • 

When decreasing Rq we find that uli decreases: the initial 
inter-particle distance gets smaller, as £ cx Rq, and thus 
hard scatterings are more probable. The dependence of 
aLi on the particle number clearly shows the effect of two- 
body collisions. For small particle number, i.e. A'' < 10^, 
the system never reaches a QSS and the LI test fails at all 
times. We then find that, at fixed e and rj, and for N £ 
[W^, 10^] the breakdown expansion factor qli is almost 
independent with N. This behaviour can be interpreted 
as due to the fact that the larger is N, the larger the 
core mass density and thus the characteristic time Tc'. in 



this situation the collisional time scale T2 in EqU] becomes 
almost independent on N. 

Finally by increasing the initial velocity dispersion, i.e. 
de creasing bg , w e find t hat qli increases: indeed, as shown 
in ISvlos Labinil (2012), when the initial configuration is 
closer to the virial equilibrium, i.e. 60 — ^Ij the collapse 
is less violent so that the size and density of the struc- 
ture after the collapse are not largely changed as it occurs 
when &o = 0. In this latter case the core of the struc- 
ture after the collapse may reach very high densities, i.e. a 
more favourable condition for the occurrence of two-body 
scatterings. However, we have found that for bo e [—1,0] 
the density profile, in the long time evolution, approaches 
the same NFW shape. Thus, the density profile, for long 
enough times, approaches the 1 /r^ behaviour for any value 
of bo a fact that shows that this tail is formed by the same 
dynamics process, independent on the initial conditions, 
in all cases. 

While the origin of the r~'^ can be understood by con- 
sidering that the energy of particles orbiting around the 
core is conserved I Svlos Labinil (j2012l) , the effect of inac- 
curately integrated two-body scatterings is indeed that 
of perturbing such energy conservation. The energy in- 
jection, is especially efficient for the fastest particles, i.e. 
those that form the tail of the density profile. A quanti- 
tative study of this effect is postponed to a forthcoming 
work, but qualitatively we may conclude that the defor- 
mation of the density profile is related to the inability 
to resolve numerically high velocity particles and thus to 
conserve their energy properly. 

7.2. Some comments on more complex initial 
conditions 

The simplicity of the system considered as initial condi- 
tions has let us to develop systematic tests to control the 
accuracy of a cosmological code. The principal result that 
we have obtained is that the range of redshifts where the 
numerical integration of a certain structure is reliable, i.e. 
it follows the stable clustering regime and it passes the LI 
test, is in general limited to narrow one that is dependent 
on the various physical (number of particles, velocity dis- 
persion, density and amplitude of the over-density) and 
numerical (softening length, time-stepping) parameters of 
a given simulation. Therefore, our results suggest that the 
limited numerical resolution of a given simulation can be 
not sufficient to integrate the formation of structures with 
a wide range in size and mass as those formed in the cur- 
rently favoured cosmological models structures through 
a hierarchical aggregation. A more detailed study of the 
more complex cosmological case is thus required. 

Here we note that, if one is not able to integrate cor- 
rectly the evolution of this simplified system, one will en- 
counter more difficult problems if one wants to properly 
simulate the formation of structures of different sizes, with 
different particle numbers and in a complex environment 
as that occurring in the case of typical initial conditions of 
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structure formation models. On the other hand the choice 
of the physical and numerical parameters adopted in this 
paper could be completely unrealistic and thus irrelevant 
for the case of a cosmological simulations that uses ini- 
tial conditions as those predicted by standard models of 
galaxy formation. 

In this respect, we note that the numerical param- 
eters have been varied in a relatively large range, i.e. 
e e [10-^,10^3] and ?7 e 2.5 X [10-3,10-1]. In addition, 
the number of particles used, ranging in N G [10'^, 10^], is 
comparable to the number of particles whi ch form a single 



halo in a typical cosmological simulation (jSpringel et al 
20051) . Thus the system we considered, once it reaches a 



virialized state, can be though to represent a simplified 
toy model approximating a virialized cosmological halo at 
its formation time. This structure, however then evolves 
as if it were isolated from the rest of the universe. Thus, 
such a self-gravitating system presents several differences 
with respect to the halos formed in typical cosmological 
simulations: it is spherical and isolated so that it has ini- 
tially a simple geometry, it is not subjected neither to mass 
accretion nor to tidal interactions with the neighbouring 
structures, it does not fragment into (macroscopic) sub- 
structures during its evolution and it does not collide with 
structures of comparable size. 

A possible method to control both the role of 
the physical and the numerical parameters in a 
more complex co smological simulations was presented in 
Jovce fc Svlos Lab ini (2013,), where it was outlined a sort 
of recipe that can be simply implemented and that can fur- 
nish a general manner for testing the accuracy of a simula- 
tion of halos containing N particles. This is simply based 
on the testing of whether an halo formed in a complex 
environment, once it is isolated, evolves in agreement with 
stable clustering and represents a quasi-stationary config- 
uration, given the same numerical parameter of the nu- 
merical integration (softening length, time-stepping, etc.) 
and the same physical parameters of the structure (num- 
ber of particles, size and velocity dispersion). The im- 
plementation of this test is postponed to a forthcoming 
work. 

As a final remark, we stress that we have observed that 
the density profile of the virialized structure formed from 
the simple initial conditions we have considered converges, 
for a wide range the numerical and physical parameters of 
the simulation, to the NFW profile, i.e. the same profile 
observed in cosmological N-body simulations that start 
from the typical initial conditions predicted by standard 
models of galaxy formation. Whether this is a coincidence 
or whether also standard cosmological halos are seriously 
affected by the numerical (in) accuracy of the numerical 
integration is an open problem that will explored in full 
details in forthcoming works 

I thank Michael Joyce for useful discussions and com- 
ments, Roberto Ammendola and Nazario Tantalo for the 
valuable assistance in the use of the Fermi supercomputer 
where the simulations have been performed. I also would 
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ful comments and suggestions that have allowed to im- 
prove the presentation of the paper. 
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Table 1. Details of the N-body simulations: N is the num- 
ber of particles in the spherical over-density, |&o| is the 
absolute value of its initial virial ratio, Rq is its initial 
radius, e is the softening length, £ is the initial average 
distance between nearest neighbours and 77 is the param- 
eter that control the accuracy of the time-stepping. Only 
in the case of the simulation M4b (M4c) we have changed 
the parameter 9 from 0.7 to 0.1 (the parameter (f) from 
5 X lO"-^ to 5 X 10^''). (See Sect 13. 21 for a discussion of the 
various numerical parameters of the code, i.e. e, rj, 6, (p). 
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5 


4/(37r)V5 


e 


0.025 


15360 


52.6 


l.OE-3 


0.05 


1920 


18.6 
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240 


6.6 


4.1E-3 


0.2 


30 


2.3 


8.2E-3 


0.3 


8.9 


1.1 


1.2E-2 



Table 2. Parameters of the spherical clouds with differ- 
ent initial radius Rq and with over-density S fEa ll2l) . The 
third column reports the pre-factor of EqHD The aver- 
age distance between nearest neighbours £ is computed 
for N = 10'*, i.e. i w 4.12 x W^'^Rq. 



